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Abstract 

We compute the pair entanglement between two interacting bosons in a two dimensional (2D) 
isotropic harmonic trap. The interaction potential is modelled by a 2D regularized pseudopoten- 
tial. By analytically decomposing the wave function into the single particle basis, we show the 
dependency of the pair entanglement on the scattering length. Our results turn out to be in good 
agreements with earlier results using a quasi-2D geometry. 
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The rapid development of laser cooling, trapping and manipulating ultracold bosonic 
atoms in optical lattices has attracted much interest due to the capability of a precise 
control over the number of atoms in each lattice site 

nasi 

|5j, thus raises hope for 

the application of optical lattice systems in quantum information and communication. Such 
lattice bosons exhibit novel properties such as superfiuid-Mott insulator quantum phase 
transition, which has been observed successfully in the laboratory On the other 

hand, recently there has bee n g rowing interest in the interacting atomic system in low 
dimensions llO[. In ultracold atomic gases, low dimensions are usually 

obtained by imposing strong confinement in some direction whose motional degree of freedom 
is frozen, thus only virtual transitions are allowed in that direction. For example, in a pioneer 
work of M. Olshanii, it has been shown that the effective one dimensional (ID) scattering 
length can be renormalized from a three dimensional (3D) scattering length opening 
intensive discussions on the subject of restricted scattering. One unique property possessed 
by ID hard core bosonic systems is the so called "Bose- Fermi duality" predicted long time 
ago jl2|, which is recently verified in ultracold atomic gases experimentally 0,0]. 

Quantum entanglement is a central concept in quantum mechanics. It is viewed as valu- 
able resources enabling the unmatched power of quantum computation and communication 
Q. In previous studies, the entanglement of .-wave scattenng between two interaeting 
)osons has been investigated in a ID harmonic trap [1^ and a 3D spherical harmonic trap 
isl ]. However, the direct calculation from the 2D pseudopotential is still absent. The recent 
work of using a 3D cylindrical harmonic trap has revealed the entanglement properties un- 
der the quasi-2D geometry jioj]. It is natural to speculate that, for very strong confinement, 
results from quasi-2D geometry calculations should give essentially the same answer as the 
true 2D case. It is the purpose of this paper to discuss the true 2D case and compare it with 
quasi-2D results. 

In this paper, we study the entanglement between two interacting bosons in a 2D isotropic 
harmonic trap. We adopt the 2D pseudopotential to model the s-wave scattering interaction. 
From the obtained spectrum, we then compute the pair entanglement for the ground state 
and the first two excited states perturbed by the interaction. Our results turn out to be 
in good agreements with earlier results using quasi-2D geometry, thus confirming the above 
speculation. 

The model we consider consists of two interacting bosons in a 2D isotropic harmonic trap. 
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The total Hamiltonian is 

2 



H = 

i=l 



where pi = (xj, r/i) is the position vector of the i-th atom in cartesian coordinate. is the 
s-wave scattering potential given by 171] 



X (^l-ln(Api2)pi2^) , (2) 

with a being the s-wave scattering length and pi2 = \pi — P2\- A is an arbitrary constant, 
which will be cancelled out if the potential acts on the proper wave functions with logarithm 
divergence (oc ln(p)). Similar to the 3D case, the regularization in the above pseudopotential 
is necessary to avoid the unwanted divergence. 

Following the usual procedure of separating the total Hamiltonian into center of mass 
(cm.) and relative (rel) motion, we write H = Hera. + -f^rei, with Hem. being a 2D har- 



monic oscillator. We further use a^o = ^/h/Jmcu) and hu as the length and energy scale, 
respectively. Thus, the corresponding coordinate transformations are g = {pi + p2)/V^ and 
p= (pi — p2)/v^, for the cm. and rel-motion respectively. The dimensionless version of 
the relative Hamiltonian is 



^rel = -^Vj.+ ^p' + K(p1, (3) 



where Vs is now given by 



For energy E = fc^, the Wigner threshold law requires fca <^ 1 in the free space. Applying 
to our case, this condition turns out to be a <^ if the threshold law still holds in an 
external trap. Otherwise, we must use an energy- dependent scattering length to obtain 
more accurate results 18|. Here we simply treat the scattering length as an independent 
parameter. Although such a simplification is supposed to be valid only for a small scattering 
length, the direct extension to a large scattering length beyond the Wigner threshold law 
allows for a closed form solution and a qualitative analysis for the current problem. While 
for an energy dependent scattering length, it can be only be obtained numerically through 
a self consistent procedure and therefore lacks an analytical treatment. 



From the relative Hamiltonian, we can see that the harmonic orbitals with angular quan- 
tum number m 7^ are obvious solutions because of the contact form of the scattering 
potential. To obtain the non-trivial solutions with scattering length dependence, we expand 
'.he eigenfunction into harmonic orbitals with vanishing angular quantum number as in Ref. 
lol . ip{p) = ^n'^n4>n{p), with the 2D s-wave harmonic orbital given by 

1 



0n(p) 



-''/'Up'). 



(5) 



Here Ln{x) is the n-th order Laguerre polynomial. The expansion coefficient takes a simple 
form as a„ oc l/(2n + 1 — E), from which we can determine the corresponding spectrum. 
After some algebra, we obtain the eigenfunction 

T(-u) 



-.e 2 



(6) 



with eigenenergy E = 21/ + 1. Here r(n) is the complete Gamma function and U{a, b, x) is 
the confluent hypergeometric function of the second kind, ({s, z) = Xl^o !/("- + is the 



Hurwitz zeta function 



20|. 



For p — > 0, we have the asymptotic results as follows [2l| 



l + 27z/ + J^(l-z/)i^ 



(7) 



r(i-^) 

where !F is the digamma function [20|. Thus the wave functions are consistent with the 
scattering potential and indeed give a result independent of A. 
The energy quantization condition is found to be 



v^exp 
-\/2exp 



-7 - 



(8) 



We note the scattering length is nonnegative [22]. The results of v as function of ln(a) are 
shown in Fig. [H We find that the energy spectrum behaves like that in a ID and 3D trap: 
(1) the energy is a smooth and monotonic function of the scattering length, and (2) the 
energy is bounded except for the ground state. 
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The radial densities plV'lp)!^ of the ground state is shown in Fig. [2] for chosen scattering 
lengths ln(a) = ±1/2 and ln(a) +00. We can see that the atoms tend to approach each 
other as In (a) decreases. In addition, we note that the derivative diverges at p = for 
ln(a) = ±1/2 which originates from the scattering potential. While for ln(a) — > ±00 it has 
a finite derivative corresponding to the non-interacting case. This is because in this limit, 
the effective coupling strength (720 — — 27r/ln(a) approaches |23l |. 

After obtaining the spectrum, we are able to compute the entanglement. We adopt von 
Neumann entropy as our pure state entanglement measure 1^. We are interested in the 
pair entanglement for pure states with cm. motion being the ground state of the trap. The 
total wave function is 



T(-u) 



(9) 



First, we briefly review the method in Ref. [l^ which is also applicable to our case. Since 

q"^ = {pj + pI + 2pip2Cos(v9i - v52))/2 and p^ = {p\ + pi - 2pip2Cos{(pi - (p2))/'2, we can 
expand the total wave function as \l>oi/(^', p) = Zlm ^m(Pi' P2)e*"''^'^^~'^^V(PiP2)- ^m(pi,P2) 
can be further written in the Schmidt form Am{pi,P2) = J2s '^m,sXs"\pi)Xs"\p2) with 

r> + 00 



AUPuP2)X^r\P2)dp2 = ^m,sx^r\pi)- 



(10) 



Thus, the total wave function can also be written in the Schmidt form 



Pi 



x^r\p.) 
( 

P2 



(11) 



The entropic entanglement can then be calculated from the Schmidt coefficients as £^ = 

Here we use an alternative method to calculate the entanglement. We find the following 
decomposition in the single particle basis 



^0i/(Pl,P2) = C'^n0m(Pl)0n(p2) 



E 



TT 



v/C(2:^t^2z.-E,(m,± 

■p f mj+Uj+l 

f-1)^^ 



<Pm,{plj)(pn,{p2j), 



(12) 



where j = {x, y}, m = {rrix, my) and n = {rix, riy). pkj is the j-th component of the position 
vector pk for the k-th atom in cartesian coordinate, e.g., pix = xi. (pmip) is the 2D harmonic 
orbital again in cartesian coordinate. This decomposition can be proved straightforwardly 
and its validity holds as long as the cm. motion is in the ground state of the trap. We 
then diagonalize the expansion coefficient matrix [C^^] and obtain U'^CU = D, where D 
is a diagonal matrix with diagonal elements A„ satisfying Yin — ^- The von Neumann 
entropy is then given by £^ = —J2n^n^^iKi)- The calculated entropy is shown in Fig. [3] 
where the (red) decreasing curve is for the ground state and the lower (blue) and upper 
(blue) increasing curves are for the first and second excited state, respectively. We find that 



the tendency of the curves remains essentially the same as in the 3D case 15|]. In addition 



the result is in a q ualitative agreement with that from a cylindrical trap with an aspect ratio 
r] = uj± <ti Uz [l(y\. It is interesting to make a quantitative comparison between a 2D and 



quasi-2D geometry. 



length is given by 17 



b do this, we note that, in the quasi-2D case, an effective scattering 
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aeff = 2.092a,exp (^-y , (13) 

where asD is the original 3D scattering length and az = ^/h/ (mujz) is the length scale along 
the tightly confined ^-axis. Transforming Eq. ( |T3l) into the unit that is adopted in the 
current paper, it reads 

In^ = In^ - ,/^^. (14) 

a_L ^ V 2 flsD 

Our comparison with 77 = 20 is also shown in Fig. [3] where data points with diamonds (cir- 
cles) are for the first (second) excited state. We do not make comparison for the ground state 



25l | . We can see that the two results agree well with each other except for the critical region 
where a^D ~ ±00 and the corresponding effective scattering length is Ina^g ■* = ln(2.092/ y/fj) 
as from Eq. f|T^ . For rj = 20, InOg^'' ^ —0.76. The 2D results vary more slowly than the 
quasi-2D results around the critical region. The discrepancy is expected since in deriving 
the effective scattering length in Eq. f|T3l) . there is no trapping potential in the transverse 
direction. While in our study, the existence of such a transverse trapping potential may 
have an effect on the scattering problem. 

To better understand the entanglement properties, we give an example for the eigenvalues 
and eigenfunctions. We randomly choose a scattering length ln(a) = —0.5359 and the 
corresponding results are shown in Figs. H] and O respectively. 



In Fig. m we show the first 20 eigenvalues for the ground state. Using a fitting curve 
A„ = an^, we obtain a ^ 0.7 and (3 ~ —0.9. This slow decrease is expected since the 
eigenfunction \I'ojy(pi, P2) is singular at the origin which can only be removed by a regularized 
pseudopotential. On the other hand, the entanglement is convergent which can be seen from 
the power law relation of the eigenvalues and further confirmed by our numerical results. 
In Fig. [5l we show the eigenf unctions (Schmidt orbitals) for four different eigenvalues for 
the ground state. The Schmidt orbital (j){p) is related to the harmonic basis by a unitary 
transformation, i.e., (j){p) = W(f){p). The p-th component of can be written as 4>^{p) = 
J2rh^i)rn^m{p)- Compared to harmonic orbitals, these orbitals experience deformation due 
to the interaction. On the other hand, the symmetry of the orbitals is still preserved due 
to the isotropic interaction. For example, if we rotate the orbital in Fig. by 90°, we 
obtain another eigenstate which is degenerate to it. These two eigenstates correspond to the 
second and third eigenvalue in Fig. HI 

Two limiting cases, i.e., ln(a) — > ±00, deserve a simple discussion. When ln(a) +00, 
we have already mentioned before that it corresponds to the non-interacting regime. While 
for ln(a) —00, we can see that the entanglement goes without limit. We identify it as 
the bound state regime, where the binding energy is calculated from Eq. ([8]) with the result 
Eb —* —2v ^ as z/ —00. The size of the bound state is on the order of a, thus 
a vanishingly small cross section in this regime. This is a quite general conclusion in the 
scattering theory. At low energy, the scattering amplitude for an /-wave is given by [26^ 

-a^^ - riky2 - i{khYk 

where 6 is a constant only depending on the scattering potential and ri is the effective range. 
The bound state can be obtained from the pole of fi. For low energy {k 0) s-wave 
scattering (/ = 0), it gives k = i/a. The corresponding wave length A ~ = a, thus a 
size of a. The bound state energy is E'^ = |A;^| = Both analysis are in agreement with 

our results. 

In conclusion, we have carried out the calculation for the entanglement between two 
interacting bosons in a 2D isotropic harmonic trap. The s-wave scattering is modelled by the 
well-known 2D regularized pseudopotential. By decomposing the wave function into single 
particle basis, we have shown the dependency of the pair entanglement on the scattering 
length. In addition, we have analyzed the Schmidt mode (natural orbit) of the system. Our 



Ji — 1 TTTZT^T' (15) 



current study only focuses on pure states, the more general entanglement properties with a 



temperature dependence will be an open question for future investigations |27|. With the 
possible realization of a 2D few-body quantum system in the future and the controllability 
of the scattering length through Feshbach resonance, we hope our results could be helpful 
to the study of quantum information and communication in low dimensional systems. 
This work is supported by grants with NSF. 



[1] Greiner M., Mandel O., Esslinger T., Hansch T. W., and Bloch I., Nature 415, (2002) 39. 
[2] Schrader D., Dotsenko I., Khudaverdyan M., Miroshnychenko Y., Rauschenbeutel A., and 

Meschede D., Phys. Rev. Lett. 93, (2004) 150501. 
[3] Maunz P., Puppe T., Schuster I., Syassen N., Pinkse P. W. H., and Rempe G., Nature 428, 

(2004) 50. 

[4] Paredes B., Widera A., Murg V., Mandel O., Foiling S., Cirac I., Shlyapnikov G., Hansch T. 

W., and Bloch I., Nature 429, (2004) 277. 
[5] Kinoshita T., Wenger T., and Weiss D., Science 305, (2004) 5687. 

[6] Kohl M., Stoeferle T., Guenter K., Koehl M., and Esslinger T., Phys. Rev. Lett. 94, (2005) 
080403; Stoeferle T., Moritz H., Guenter K., Koehl M., and Esslinger T., Phys. Rev. Lett. 96, 

(2006) 030401. 

[7] Sheehy D. E. and Radzihovsky L., Phys. Rev. Lett. 95, (2005) 130401. 
[8] Orignac E. and Citro R., Phys. Rev. A 73, (2006) 063611. 

[9] Petrov D.S., Holzmann M., and Shlyapnikov G. V., Phys. Rev. Lett 84, (2000) 2551; Fedichev 
P. O., Bijlsma M. J., and Zoller P., Phys. Rev. Lett. 92, (2004) 080401. 
[10] Astrakharchik G. E., Boronat J., Kurbakov I. L., and Lozovik Yu. E., Phys. Rev. Lett. 98, 

(2007) 060405. 

[11] Olshanii M., Phys. Rev. Lett. 81, (1998) 938. 

[12] Girardeau M., J. Math. Phys. (N.Y.) 1, (1960) 516. 

[13] Nielsen M. A. and Chuang I. S., Quantum computation and quantum information, Cambridge 

University Press (2000). 
[14] Sun B., Zhou D. L., and You L., Phys. Rev. A 73, (2006) 012336. 
[15] Wang J., Law C. K., and Chu M.-C., Phys. Rev. A 72, (2005) 022346. 



8 



[16] Sun B., You L., and Zhou D. L., Phys. Rev. A 75, (2007) 012332. 

[17] Pricoupenko L. and Olshanii M., J. Phys. B: At. Mol. Opt. Phys. 40, (2007) 2065. 

[18] Bolda E. L., Tiesinga E., and Julienne P. S., Phys. Rev. A 66, (2002) 013403. 

[19] Busch T., Englert B.-G., Rzazewski K., and Wilkens M., Found. Phys. 28, (1998) 549. 

[20] Arfken G., Mathematical Methods for Physicists, 3rd Ed., Academic Press, Inc., (1987). 

[21] For p^O, we can use dC/(-i^,l,p2)/dp2 = i,u{-u + 1,2, p^) u[l/(T{l - u)p'^) + 0(1)] and 
integrate on both sides to obtain U{—i', 1, p^) — > Bi,\n.[p) + Ci,- Then can be obtained by us- 
ing the Wronskian of two confluent hypergeometric functions M{—v,\, p'^)dU{—v, 1, fp'^jdfp' — 
UdM/dp^ —(1 + l//9^)/r(— z^). Comparing both sides, we immediately obtain Eq. ([7]). 

[22] Verhaar B. J., Eijnde J. P. H. W. van den, Voermans M. A. J., and Schaffrath M. M. J., J. 
Phys. A 17, (1984) 595. 

[23] Idziaszek Z. and Calarco T., Phys. Rev. A 74, (2006) 022712. 

[24] Petrov D. S. and Shlyapnikov G. V., Phys. Rev. A 64, (2001) 012706. 

[25] The ground state entanglement in a quasi-2D geometry is computationly intensive. Even for a 
moderate scattering length, it requires much more basis for the expansion than the 2D case in 
order to achieve convergence. Due to this numerical difficulty, the ground state entanglement 
is not given in Ref. [l^ and as a result, a comparison is not carried out. 

[26] Landau L. D. and Lifschitz E. M., Vol. 3 Quantum mechanics, Nonrelativistic theory (3ed., 
Pergamon, 1991). 

[27] For the thermal entanglement which is related to mixed states in general, we can use negativ- 



ity as entanglement measure 



281 ]. Different from pure state calculation which only need the 



reduced density matrix, the calculation of negativity need the full two-body density matrix. 
This is a computationally demanding task due to the very large Hilbert space. 
[28] Vidal G. and Werner R. F., Phys. Rev. A 65, (2002) 032314. 



9 



-1 

-2 
-3 



ln(a) 



FIG. 1: Energy spectrum v as function of ln(a). The unmarked red curves are for those solutions 
affected by the scattering potential, and the marked blue horizontal lines are for asymptotic values 
as ln(a) ±cxd. 
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FIG. 2: The radial density distribution of the ground state. 
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FIG. 3: Von Neumann entropy as function of the scattering length. The (red) decreasing curve 
is for the ground state. The lower (blue) and upper (black) of the two increasing curves are for 
the first and second excited state (both perturbed by interaction) in the relative motion. The cm. 
motion is always assumed in the ground state of the trap. Data points with circles (blue) and 
diamonds (black) are from quasi-2D calculations with rj = 20 and an effective scattering length of 
Eq. ([13]). See text for details. 
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FIG. 4: The first 20 eigenvalues A„ in descending order for the ground state with ln(a) = —0.5359. 
The (red) dashed dotted line is an exponential fit. 
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FIG. 5: Several Schmidt orbitals for the ground state with ln(a) = —0.5359. From (a)^(d), the 
corresponding eigenvalues are 0.7408,0.3164,0.1686, and 0.1029, respectively. 
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